Sandbox/Old plotit/plot_trafolm.R

In akreutzmann/trafo: Estimation, Comparison and Selection of Transformations

#' Plot for optimal transformation parameter - linear models
#' 
#' @param lambdarange a numeric vector with two elements defining an interval 
#' that is used for the estimation of the optimal transformation parameter.
#' @param lambdaoptim optimal or given transformation parameter.
#' @param measoptim measure at the optimal transformation parameter.
#' @param y dependent variable.
#' @param x matrix of regressors
#' @param method a character string. In order to determine the optimal parameter
#' for the transformation five different estimation methods can be chosen
#' (i) Maximum-Likelihood ("ml"); (ii) skewness minimization ("skew");
#' (iii) minimization of Kolmogorov-Smirnoff divergence  ("div.ks");
#' (iv) minimization of Craemer von Mises divergence ("div.cvm"); (v)
#' minimization of Kullback Leibner divergence  ("div.kl"). In case of no and
#' log transformation "NA" can be selected since no optimization is neccessary
#' for these two transformation types.
#' @param trafo a character string that selects the transformation.
#' @keywords internal


plot_trafolm <- function(lambdarange, lambdaoptim, measoptim, 
                         y, x, trafo, method) {
  
  lambdavector <- seq(lambdarange[1], lambdarange[2], 0.025)
  l <- length(lambdavector)
  lambdavector[l + 1]  <- lambdaoptim
  lambdavector <- sort(lambdavector)
  measvector <- sapply(lambdavector, estim_lm, y = y, x = x, trafo = trafo,
                       method = method)
  vline <- lambdaoptim
  
  if (method == "ml" | method == "reml") {
    measvector <- -measvector
    data1 <- data.frame(measvector = measvector,  lambdavector = lambdavector)  
    measoptim <- -measoptim
    y_lab <- "Profile log-likelihood"
    
    
  } else if (method == "skew" | method == "pskew") {
    data1 <- data.frame(measvector = measvector,  lambdavector = lambdavector)
    y_lab <- "Skewness"
  } else if (method == "div.ks" | method == "div.cvm" | method == "div.kl") {
    data1 <- data.frame(measvector = measvector,  lambdavector = lambdavector)
    y_lab <- "Divergence"
  }
  
  #plot <- ggplot(data1, aes(x = lambdavector,
  #                          y = measvector)) + geom_line() + 
  #  geom_vline(xintercept = vline, linetype = "dashed") + 
  #  geom_hline(yintercept = measoptim, color = "red", linetype = "dashed") + 
  #  xlab(expression(lambda)) + ylab(y_lab)
  
  dev.hold()
  plot(data1$lambdavector, data1$measvector, type = "l", lwd = 1.5,
      xlab = expression(lambda), ylab = y_lab) 
  abline(h = measoptim, lty = 2, col = "red")
  abline(v = lambdaoptim, lty = 2)
  dev.flush()


  out <- list(lambdavector = lambdavector, 
              measvector = measvector)
  return(out)
  }